Time Series Final Group Project
Mia Song
2023-05-22
Import necessary packages
library(TSA)##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(fpp)## Loading required package: forecast## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.Arima TSA
## plot.Arima TSA## Loading required package: fma## Loading required package: expsmooth## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: tserieslibrary(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(readxl)
library(corrplot)## corrplot 0.92 loadedlibrary(stats)1. Explanatory Analysis
1.1. Import data
df_raw <- read_excel('TS_Project.xlsx')
df <- df_raw1.2. Check data structure and nulls
head(df)## # A tibble: 6 × 9
## time export import cpi ppi bond sentiment USDX uncertai…¹
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2002-01-01 00:00:00 1.41 3.57 178. 128. 5.04 93 120. 117.
## 2 2002-02-01 00:00:00 1.34 3.52 178 128. 4.91 90.7 119. 87.6
## 3 2002-03-01 00:00:00 1.76 3.79 178. 130. 5.28 95.7 119. 83.4
## 4 2002-04-01 00:00:00 1.46 4.45 179. 131. 5.21 93 115. 88.6
## 5 2002-05-01 00:00:00 1.63 4.36 180. 131. 5.16 96.9 112. 86.3
## 6 2002-06-01 00:00:00 1.76 4.48 180. 131. 4.93 92.4 106. 93.6
## # … with abbreviated variable name ¹​uncertaintytail(df)## # A tibble: 6 × 9
## time export import cpi ppi bond sentiment USDX uncertai…¹
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2022-09-01 00:00:00 6.72 28.5 297. 268. 3.52 58.6 112. 174.
## 2 2022-10-01 00:00:00 5.91 28.4 298. 265. 3.98 59.9 112. 177.
## 3 2022-11-01 00:00:00 5.50 25.6 299. 263. 3.89 56.8 106. 172.
## 4 2022-12-01 00:00:00 5.85 24.8 299. 258. 3.62 59.7 104. 136.
## 5 2023-01-01 00:00:00 6.16 25.0 301. 260. 3.53 64.9 102. 143.
## 6 2023-02-01 00:00:00 5.86 20.6 302. 259. 3.75 67 105. 124.
## # … with abbreviated variable name ¹​uncertaintyany(is.na(df))## [1] FALSE1.3. Transform as Time Series data
export = ts(df$export, start = c(2002, 1), frequency = 12)
import = ts(df$import, start = c(2002, 1), frequency = 12)
cpi = ts(df$cpi, start=c(2002,1),frequency=12)
ppi = ts(df$ppi, start=c(2002,1),frequency=12)
bond = ts(df$bond, start=c(2002,1),frequency=12)
sentiment = ts(df$sentiment, start=c(2002,1),frequency=12)
usdx = ts(df$USDX, start=c(2002,1),frequency=12)
uncertainty = ts(df$uncertainty, start=c(2002,1),frequency=12)1.4. Correlation
# create a correlation matrix
cor_matrix <- cor(cbind(export, import,
cpi, ppi,
bond, sentiment,
usdx, uncertainty))
# Plot a correlation heatmap
corrplot(cor_matrix, method = "color",
order = "hclust",
tl.col = "black", tl.srt = 45, addCoef.col = "black")
# Print the correlation table
# CPI & PPI highly correlated, only keep one.
print(cor_matrix)## export import cpi ppi bond
## export 1.0000000 0.9576931 0.9508205 0.88258986 -0.61290296
## import 0.9576931 1.0000000 0.9623065 0.91101757 -0.64355565
## cpi 0.9508205 0.9623065 1.0000000 0.92833707 -0.69498684
## ppi 0.8825899 0.9110176 0.9283371 1.00000000 -0.56578495
## bond -0.6129030 -0.6435556 -0.6949868 -0.56578495 1.00000000
## sentiment -0.1779160 -0.2435011 -0.2186952 -0.43137804 0.07714775
## usdx 0.2673713 0.2693136 0.2465156 0.01599541 -0.07202055
## uncertainty 0.3651635 0.4453334 0.4491599 0.38277243 -0.59587323
## sentiment usdx uncertainty
## export -0.17791601 0.26737131 0.36516347
## import -0.24350112 0.26931356 0.44533338
## cpi -0.21869519 0.24651557 0.44915994
## ppi -0.43137804 0.01599541 0.38277243
## bond 0.07714775 -0.07202055 -0.59587323
## sentiment 1.00000000 0.32526534 -0.46620569
## usdx 0.32526534 1.00000000 0.04882897
## uncertainty -0.46620569 0.04882897 1.000000001.5. Time Series plots for each variable
autoplot(import, main="Import") + ylab('Import(in billion $)')
autoplot(cpi, main="CPI") + ylab('CPI')
autoplot(ppi, main="PPI") + ylab('PPI')
autoplot(bond, main="10Y Govn't Bond") + ylab('Bond')
autoplot(sentiment, main='Consumer Sentiment') + ylab('Sentiment')
autoplot(usdx, main='USD Index') + ylab('USD Index')
autoplot(uncertainty, main='Economic Uncertainty Index') + ylab('Economic Uncertainty Index')
1.6. Stationarity
1.6.1. BoxCox transformation
lambda <- BoxCox.lambda(import)
lambda## [1] -0.08103269cbind('Before Transformation' = import,
'After Transformation'= BoxCox(import, lambda = lambda)) %>% autoplot(facet = TRUE) + ggtitle('BoxCox Transformation Comparison') + ylab('Import(B$)')
1.6.2. Differencing & KPSS Test
# KPSS Test with the original data
kpss.test(import)## Warning in kpss.test(import): p-value smaller than printed p-value##
## KPSS Test for Level Stationarity
##
## data: import
## KPSS Level = 3.7518, Truncation lag parameter = 5, p-value = 0.01# KPSS Test with 1st order differencing data
diff1 = diff(import)
kpss.test(diff1)## Warning in kpss.test(diff1): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: diff1
## KPSS Level = 0.027579, Truncation lag parameter = 5, p-value = 0.1plot(diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing")
# KPSS Test with 1st order differncing + baxcox data
bc_diff1 = diff(BoxCox(import, lambda = lambda))
kpss.test(bc_diff1)## Warning in kpss.test(bc_diff1): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: bc_diff1
## KPSS Level = 0.10007, Truncation lag parameter = 5, p-value = 0.1plot(bc_diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing with BoxCox Transformation")
1.6.3. Differencing & ACF, PACF Plot
# Original Data
Acf(import)
Pacf(import)
# Data with 1st order differencing
Acf(diff1)
Pacf(diff1)
# Data with 1st order differencig + boxcox
Acf(bc_diff1)
Pacf(bc_diff1)
1.7. Check seasonlaity
- We perform seasonlaity insepctio after making our data stationary, to remove other components to spot seasonlaity in our data more clearly.
# Original Plot
autoplot(import)
# Decomposition Plot (STL)
ts_data <- ts(bc_diff1, frequency = 12)
stl_fit <- stl(ts_data, s.window = "periodic")
plot(stl_fit)
remainder <- stl_fit$time.series[, "remainder"]
plot(remainder)
# Seasonal Plot
ggseasonplot(bc_diff1, polar = TRUE)
ggseasonplot(bc_diff1, polar = FALSE)
# Subseries Plot
ggsubseriesplot(bc_diff1) +
labs(y = "Import(B$)",
title = "Seasonal plot: Import")
# Periodogram
spectrum(bc_diff1)
periodogram(bc_diff1)
temp <- periodogram(bc_diff1)
max_freq <- temp$freq[which.max(temp$spec)]
seasonality <- 1/max_freq
1 / temp$freq[43] # 5.953488## [1] 5.9534881 / temp$freq[107] #2.392523## [1] 2.3925231 / temp$freq[21] # 12.19048## [1] 12.19048Acf(bc_diff1)
Pacf(bc_diff1)
2. Model Exploration: ETS
- From this part, we will redefine our df, leaving only ‘import’ feature.
- We will first fit our model using the period from Jan/2002 ~ Dec/2018 to predict 2019. Then, we will gain fit our model using the period from Jan/2002 ~ Dec/2021 to predict 2022. This is to compare the impact of COVID19 for our time series analysis and the best model.
Before COVID 19
2.1. Data preprocessing
2.1.1. Import data
df_raw <- read_excel('TS_Project.xlsx')
df <- ts(df_raw$import, start=c(2002,01), frequency=12)2.1.2. Plot data
autoplot(df) + ylab('Import') + ggtitle('US Monthly Import(2002/01 ~ 2023/02)')
2.1.3. Train-Test Split
df_train <- window(df, start=c(2002, 1), end=c(2018, 12))
df_test <- window(df, start=c(2019, 1), end=c(2019, 12))2.2. ETS Modeling
2.2.0 Seasonal Naive(Baseline Model)
fcst_snaive <- snaive(df_train, h = 12)
autoplot(snaive(df_train, h = 12))
accuracy(fcst_snaive)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.7617909 1.527044 1.280903 7.1367 12.59435 1 0.7750122accuracy(fcst_snaive$mean, df_test)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.4499564 1.248044 1.03557 2.415875 5.877225 -0.3960617 0.7297898checkresiduals(fcst_snaive) 
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 682.64, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 242.2.1 Simple Exponential Smoothing
- Simple exponential smoothing is a time series forecasting method that uses a weighted average of past observations to make future predictions.
- The method is based on the assumption that the future values of a time series can be estimated as a weighted average of the past values, with more recent observations being given greater weight.
ses3 <- ets(df_train, model='ANN')
fcst_ses3 <- forecast(ses3, h=12)
accuracy(fcst_ses3)[,2] ## [1] 0.9509481accuracy(fcst_ses3$mean, df_test)[,2] ## [1] 1.643461checkresiduals(fcst_ses3$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 172.5, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24autoplot(ses3) + ylab('Import($B)')
autoplot(ses3) + autolayer(fitted(ses3)) + ylab('Import($B)')
2.2.2. Holt’s Exponential Smoothing
- The ETS model assumes that the observed time series can be decomposed into three components: trend, seasonal, and error.
- The trend component refers to the long-term pattern or direction of the time series. The seasonal component refers to the repeating pattern that occurs at fixed intervals of time, such as daily, weekly, or yearly. The error component refers to the random and unpredictable variation in the data that cannot be explained by the trend or seasonal components.
# 1) Linear trend with additive seasonality
hwa <- hw(df_train, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
## ACF1
## Training set 0.01406249accuracy(hwa$mean, df_test) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.1271767 0.8469332 0.6890742 0.4701921 3.986005 -0.4069091 0.4849101checkresiduals(hwa$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
##
## Model df: 0. Total lags used: 24# pdf("plot.pdf", width = 10, height = 6)
plot(hwa)
lines(df_test, col = "red")
# 2) Linear trend with multiplicative seasonality
hwm <- hw(df_train, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.02226477 0.5818108 0.4363622 -0.01005463 4.344781 0.3406676
## ACF1
## Training set -0.03132167accuracy(hwm$mean, df_test) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.2187908 0.8660047 0.7331691 0.9624011 4.256856 -0.3067592 0.5090764checkresiduals(hwm$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 42.431, df = 24, p-value = 0.01154
##
## Model df: 0. Total lags used: 24plot(hwm)
lines(df_test, col = "red")
# 3) Linear trend with additive seasonality and damping
hwad <- hw(df_train, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.08943375 0.6597823 0.4971829 0.8109482 5.319108 0.3881503
## ACF1
## Training set -0.01316333accuracy(hwad$mean, df_test) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.2356866 1.105793 0.913462 0.7802561 5.422637 0.1745505 0.6544695checkresiduals(hwad$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.299, df = 24, p-value = 0.01986
##
## Model df: 0. Total lags used: 24plot(hwad)
lines(df_test, col = "red")
# 4) Linear trend with multiplicative seasonality and damping
hwmd <- hw(df_train, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.03776369 0.5807925 0.4377261 0.2958402 4.469252 0.3417324
## ACF1
## Training set -0.02054866accuracy(hwmd$mean, df_test) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.4950831 0.9971667 0.8049686 2.566052 4.561863 -0.2181502 0.5892602checkresiduals(hwmd$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 35.775, df = 24, p-value = 0.0577
##
## Model df: 0. Total lags used: 24plot(hwmd)
lines(df_test, col = "red")
2.2.3. Holt-Winter’s Exponenetial Smoothing
# 1 ETS_BEST
ets_best <- ets(df_train)
summary(ets_best) ## ETS(M,A,M)
##
## Call:
## ets(y = df_train)
##
## Smoothing parameters:
## alpha = 0.5353
## beta = 1e-04
## gamma = 2e-04
##
## Initial states:
## l = 4.0986
## b = 0.0424
## s = 0.9908 1.0682 1.1106 1.0167 1.0183 1.0399
## 1.0286 0.9756 0.9869 0.9801 0.8526 0.9316
##
## sigma: 0.057
##
## AIC AICc BIC
## 847.9923 851.2826 904.4003
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03831786 0.5906786 0.4378753 0.1381647 4.319084 0.3418489
## ACF1
## Training set -0.03521243checkresiduals(ets_best$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 42.295, df = 24, p-value = 0.01196
##
## Model df: 0. Total lags used: 24fcst_ets_best <- forecast(ets_best, h=12)
accuracy(fcst_ets_best$mean, df_test)[,2] #RMSE:2.534726## [1] 0.8525379autoplot(forecast(ets_best)) + ylab('Import($B)')
plot(ets_best)
plot(fcst_ets_best)
lines(df_test, col = "red")
# 2. ETS_AAA
ets1 <- ets(df_train, model='AAA')
summary(ets1) ## ETS(A,A,A)
##
## Call:
## ets(y = df_train, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.4797
## beta = 1e-04
## gamma = 0.2997
##
## Initial states:
## l = 4.0752
## b = 0.0591
## s = -0.2185 0.1891 0.3965 -0.0095 0.2708 0.3481
## 0.3534 0.1724 0.1472 -0.1551 -0.8346 -0.6597
##
## sigma: 0.6759
##
## AIC AICc BIC
## 942.3971 945.6874 998.8051
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
## ACF1
## Training set 0.01406249checkresiduals(ets1$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
##
## Model df: 0. Total lags used: 24fcst_ets1 <- forecast(ets1, h=12)
accuracy(fcst_ets1$mean, df_test)[,2] ## [1] 0.8469332autoplot(forecast(ets1)) + ylab('Import($B)')
plot(ets1)
plot(fcst_ets1)
lines(df_test, col = "red")
# 3. ETS MAA
ets2 <- ets(df_train, model='MAA')
summary(ets2) ## ETS(M,A,A)
##
## Call:
## ets(y = df_train, model = "MAA")
##
## Smoothing parameters:
## alpha = 0.5215
## beta = 1e-04
## gamma = 0.206
##
## Initial states:
## l = 4.1148
## b = 0.0575
## s = 0.0136 0.3469 0.5217 0.1403 0.2663 0.2594
## 0.1876 -0.0763 -0.0638 -0.1533 -0.8209 -0.6217
##
## sigma: 0.0615
##
## AIC AICc BIC
## 880.9684 884.2588 937.3765
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.007855664 0.6566577 0.4792403 -0.2734595 4.698261 0.3741425
## ACF1
## Training set 0.004567237checkresiduals(ets2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 38.875, df = 24, p-value = 0.02815
##
## Model df: 0. Total lags used: 24fcst_ets2 <- forecast(ets2, h=12)
accuracy(fcst_ets2$mean, df_test)[,2] ## [1] 0.8176148autoplot(forecast(ets2)) + ylab('Import($B)')
plot(ets2)
plot(fcst_ets2)
lines(df_test, col = "red")
# 4. ETS MMM
ets3 <- ets(df_train, model='MMM')
summary(ets3) ## ETS(M,M,M)
##
## Call:
## ets(y = df_train, model = "MMM")
##
## Smoothing parameters:
## alpha = 0.554
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 4.037
## b = 1.0072
## s = 0.9947 1.0738 1.1076 1.0177 1.0197 1.0366
## 1.0237 0.9757 0.9852 0.9811 0.8512 0.933
##
## sigma: 0.0567
##
## AIC AICc BIC
## 847.6910 850.9813 904.0990
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01799567 0.5900927 0.4405025 -0.2720982 4.360535 0.3438999
## ACF1
## Training set -0.05031675checkresiduals(ets3) 
##
## Ljung-Box test
##
## data: Residuals from ETS(M,M,M)
## Q* = 40.548, df = 24, p-value = 0.01866
##
## Model df: 0. Total lags used: 24fcst_ets3 <- forecast(ets3, h=12)
accuracy(fcst_ets3$mean, df_test)[,2] ## [1] 0.8680253autoplot(forecast(ets3)) + ylab('Import($B)')
plot(ets3)
plot(fcst_ets3)
lines(df_test, col = "red")
After COVID19
3.1. Data preprocessing
3.1.1. Train-Test Split
df_train2 <- window(df, start=c(2002, 1), end=c(2021, 12))
df_test2 <- window(df, start=c(2022, 1), end=c(2022,12))3.1.2 Seasonal Naive(Baseline Model)
fcst_snaive <- snaive(df_train2, h = 12)
autoplot(snaive(df_train2, h = 12))
accuracy(fcst_snaive $mean, df_test2)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 3.242702 4.445163 4.111537 12.07819 15.53626 0.4679377 1.6112653.2. ETS Modeling
3.2.1 Simple Exponential Smoothing
ses32 <- ets(df_train2, model='ANN')
fcst_ses32 <- forecast(ses32, h=12)
accuracy(fcst_ses32)[,2] ## [1] 1.118179accuracy(fcst_ses32$mean, df_test2)[,2] ## [1] 2.103524checkresiduals(fcst_ses32$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 211.62, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24autoplot(ses32) + ylab('Import($B)')
autoplot(ses32) + autolayer(fitted(ses32)) + ylab('Import($B)')
3.2.2. Holt’s Exponential Smoothing
# 1) Linear trend with additive seasonality
hwa2 <- hw(df_train2, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa2) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
## ACF1
## Training set 0.01213894accuracy(hwa2$mean, df_test2) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.142104 3.432513 2.447793 -8.391731 9.549002 0.4576512 1.273803checkresiduals(hwa2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
##
## Model df: 0. Total lags used: 24plot(hwa2)
lines(df_test2, col = "red")
# 2) Linear trend with multiplicative seasonality
hwm2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm2)## ME RMSE MAE MPE MAPE MASE
## Training set 0.04448036 0.7014787 0.5107636 0.005129007 4.569962 0.3460351
## ACF1
## Training set 0.06717551accuracy(hwm2$mean, df_test2) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.570837 3.174066 2.495035 -6.166118 9.709162 0.4946886 1.172481checkresiduals(hwm2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 60.478, df = 24, p-value = 5.467e-05
##
## Model df: 0. Total lags used: 24plot(hwm2)
lines(df_test2, col = "red")
# 3) Linear trend with additive seasonality and damping
hwad2 <- hw(df_train2, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad2) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.07959986 0.7972801 0.6008245 0.5927588 5.847758 0.4070501
## ACF1
## Training set 0.05177696accuracy(hwad2$mean, df_test2) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -3.092347 3.829269 3.108496 -12.07427 12.13198 0.2978517 1.418878checkresiduals(hwad2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 58.148, df = 24, p-value = 0.0001161
##
## Model df: 0. Total lags used: 24plot(hwad2)
lines(df_test2, col = "red")
# 4) Linear trend with multiplicative seasonality and damping
hwmd2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd2) ## ME RMSE MAE MPE MAPE MASE
## Training set 0.05315586 0.6975741 0.5102869 0.1244391 4.548272 0.3457122
## ACF1
## Training set 0.02301767accuracy(hwmd2$mean, df_test2) ## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.176211 3.633945 2.760502 -8.450697 10.74315 0.5221199 1.339611checkresiduals(hwmd2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 52.418, df = 24, p-value = 0.0006899
##
## Model df: 0. Total lags used: 24plot(hwmd2)
lines(df_test2, col = "red")
3.2.3. Holt-Winter’s Exponenetial Smoothing
# 1. ETS_BEST
ets_best2 <- ets(df_train2)
summary(ets_best2) ## ETS(M,A,M)
##
## Call:
## ets(y = df_train2)
##
## Smoothing parameters:
## alpha = 0.5311
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 4.0927
## b = 0.0597
## s = 1.0063 1.0846 1.1101 1.0194 1.0135 1.0378
## 1.0187 0.9762 0.9892 0.9744 0.8362 0.9336
##
## sigma: 0.0581
##
## AIC AICc BIC
## 1094.530 1097.287 1153.701
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06211406 0.7294272 0.5263934 0.0311818 4.498561 0.3566241
## ACF1
## Training set 0.01879104checkresiduals(ets_best2$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 46.18, df = 24, p-value = 0.004212
##
## Model df: 0. Total lags used: 24fcst_ets_best2 <- forecast(ets_best2, h=12)
accuracy(fcst_ets_best2$mean, df_test2)[,2] ## [1] 2.019115autoplot(forecast(ets_best2)) + ylab('Import($B)')
plot(ets_best2)
plot(fcst_ets_best2)
lines(df_test2, col = "red")
# 2. ETS_AAA
ets12 <- ets(df_train2, model='AAA')
summary(ets12) ## ETS(A,A,A)
##
## Call:
## ets(y = df_train2, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.4867
## beta = 0.0276
## gamma = 0.346
##
## Initial states:
## l = 4.2205
## b = 0.1176
## s = 0.202 0.0634 0.5544 -0.1505 -0.2504 0.3038
## 0.2454 0.0457 0.3895 0.1847 -0.6942 -0.8938
##
## sigma: 0.7971
##
## AIC AICc BIC
## 1223.920 1226.677 1283.091
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
## ACF1
## Training set 0.01213894checkresiduals(ets12$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
##
## Model df: 0. Total lags used: 24fcst_ets12 <- forecast(ets12, h=12)
accuracy(fcst_ets12$mean, df_test2)[,2] ## [1] 3.432513autoplot(forecast(ets12)) + ylab('Import($B)')
plot(ets12)
plot(fcst_ets12)
lines(df_test2, col = "red")
# 3. ETS MAA
ets22 <- ets(df_train2, model='MAA')
summary(ets22) ## ETS(M,A,A)
##
## Call:
## ets(y = df_train2, model = "MAA")
##
## Smoothing parameters:
## alpha = 0.5565
## beta = 1e-04
## gamma = 0.2409
##
## Initial states:
## l = 4.3895
## b = 0.0618
## s = -0.042 0.2136 0.6139 0.1723 0.2651 0.2256
## 0.2715 -0.1323 0.0125 -0.2621 -0.7543 -0.5837
##
## sigma: 0.062
##
## AIC AICc BIC
## 1126.354 1129.110 1185.524
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05355054 0.773105 0.5608316 -0.1563987 4.828425 0.3799555
## ACF1
## Training set 0.01192749checkresiduals(ets22$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 35.481, df = 24, p-value = 0.06159
##
## Model df: 0. Total lags used: 24fcst_ets22 <- forecast(ets22, h=12)
accuracy(fcst_ets22$mean, df_test2)[,2] ## [1] 1.907392autoplot(forecast(ets22)) + ylab('Import($B)')
plot(ets22)
plot(fcst_ets22)
lines(df_test2, col = "red")
# 4. ETS MMM
ets32 <- ets(df_train2, model='MMM')
summary(ets32) ## ETS(M,M,M)
##
## Call:
## ets(y = df_train2, model = "MMM")
##
## Smoothing parameters:
## alpha = 0.5646
## beta = 1e-04
## gamma = 9e-04
##
## Initial states:
## l = 4.0127
## b = 1.0074
## s = 1.0144 1.0841 1.1101 1.0237 1.0135 1.0326
## 1.0172 0.9665 0.9824 0.9872 0.838 0.9304
##
## sigma: 0.058
##
## AIC AICc BIC
## 1094.544 1097.301 1153.715
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01542948 0.7196395 0.5242459 -0.1484539 4.482695 0.3551692
## ACF1
## Training set -0.03362053checkresiduals(ets32$residuals) 
##
## Ljung-Box test
##
## data: Residuals
## Q* = 46.808, df = 24, p-value = 0.003535
##
## Model df: 0. Total lags used: 24fcst_ets32 <- forecast(ets32, h=12)
accuracy(fcst_ets32$mean, df_test2)[,2] ## [1] 2.680248autoplot(forecast(ets32)) + ylab('Import($B)')
plot(ets32)
plot(fcst_ets32)
lines(df_test2, col = "red")